! Constant and State Equations for Perfect Gas

module DataPerfectGas

    real(8) :: Gamma = 1.4d0
    real(8) :: cv = 717.65d0
    real(8) :: cp = 1004.71 ! Gamma * cv
    real(8) :: nu
    real(8) :: BulkViscConst = 1d0 ! BulkVisc = BulkViscCons * nu
    real(8) :: BulkVisc = 1.5d-5 ! bulk viskocity
    real(8) :: BulkViscTotal ! BulkViscTotal = BulkVisc - 2/3 nu
    
    
    real(8) :: GammX     ! (Gamma-1)/2Gamma
    real(8) :: invGammX     ! 2Gamma/(Gamma-1)
    real(8) :: GammY     ! 2sqrt(Gamma)/(Gamma-1)
    real(8) :: GammPower ! 1/2Gamma
    real(8) :: invGamma
    
    ! Dynamic viscosity defined by Sutherland's Law
    ! Visc = Visc(T)
    ! Visc = ViscRef * (T/TempRef)**1.5 * (TempRef+SconstVisc) / (T+SconstVisc)
    logical :: SutherlandsLaw = .false. ! flag - to use Sutherland's Law
    real(8) :: ViscRef = 1.72d-05 ! kg/(ms)
    real(8) :: TempRef = 273.15 ! K
    real(8) :: SconstVisc = 110.4 ! constant for Sutherland's Law

contains

    subroutine CalcGammas()
    implicit none
    
        invGamma = 1d0 / (Gamma)                ! 1/Gamma
        GammPower = 0.5d0 / Gamma               ! 1/2Gamma
        GammX = (Gamma-1d0) * 0.5d0 / Gamma     ! (Gamma-1)/2Gamma
        invGammX = 1d0 / GammX                  ! 2Gamma/(Gamma-1)
        GammY = 2d0 * sqrt(Gamma) / (Gamma-1d0)
    
    end subroutine CalcGammas
    
    real(8) function EntropyPRho(Pres,Rho)
    implicit none
        real(8),intent(in) :: Pres, Rho
    
        EntropyPRho = dlog(Pres) - Gamma * dlog(Rho)
    
    end function EntropyPRho

    real(8) function SoundVelocityPRho(Pres,Rho)
    implicit none
        real(8),intent(in) :: Pres, Rho
        
        SoundVelocityPRho = dsqrt(Gamma / Rho * Pres)

    end function SoundVelocityPRho


    real(8) function TemperaturePRho(Pres,Rho)
    implicit none
        real(8),intent(in) :: Pres, Rho
        
        TemperaturePRho = Pres * (0.4d0 * cv * Rho)**(-1)
        
    end function TemperaturePRho

    real(8) function PressureAllVars(Vars)
    implicit none
        real(8),intent(in) :: Vars(5) ! 1:5 = Rho, Velocities, Energy
        
        real(8) :: Rho,U(3),E
        
        Rho = Vars(1)
        U = Vars(2:4)
        E = Vars(5)

        PressureAllVars =  (E-0.5d0*sum(U(1:3)*U(1:3))) * Rho * (Gamma-1.0d0)
!        PressureAllVars = (E-0.5d0*sum(U(1:3)*U(1:3))) * Rho * 0.4d0

    end function PressureAllVars

    real(8) function DensityPS(Pres,Entropy)
    implicit none
        real(8),intent(in) :: Pres,Entropy

        DensityPS = (Pres/exp(Entropy))**invGamma

    end function DensityPS
    
    real(8) function DensityPT(Pres,Temp)
    implicit none
        real(8),intent(in) :: Pres,Temp
        
        DensityPT = Pres * (Temp*Cv*(Gamma-1.0d0))**(-1)

    end function DensityPT

    real(8) function EnergyPRhoU(Pres,Rho,U)
    implicit none
        real(8),intent(in) :: Pres,Rho,U(3)
        
        EnergyPRhoU = Pres * ((Gamma-1.0d0)*Rho)**(-1) + 0.5d0 * sum(U(1:3)*U(1:3))

    end function EnergyPRhoU

    real(8) function ACoefPRho(Pres,Rho)
    implicit none
        real(8),intent(in) :: Pres,Rho

        real(8) :: s

!        ACoefPRho = GammY * (Pres * Rho**(-Gamma))**GammPower
        
        s = dlog(Pres)-Gamma*dlog(Rho)
        ACoefPRho = GammY * dexp(s*GammPower)
        
    end function ACoefPRho
    
    real(8) function DynamicViscosity(Temp, Rho)
    implicit none
        real(8),intent(in) :: Temp, Rho

        if(SutherlandsLaw)then
            DynamicViscosity = ViscRef * (Temp/TempRef)**1.5 * (TempRef+SconstVisc) / (Temp+SconstVisc)
        else
            DynamicViscosity = Rho * nu ! nu !
        end if

    end function DynamicViscosity
    
    real(8) function TotalBulkViscosity(DynVisc)
    use Data_cfoam, only: ZeroEpsilon
    implicit none
        real(8), intent(in) :: DynVisc
    
        if(DynVisc.gt.ZeroEpsilon)then
	        TotalBulkViscosity = DynVisc * (BulkViscConst - 2d0 / 3d0) * 0.5
	                ! multiplicator 0.5 because bulk viscosity added to Grad_i U_j
	                ! where gradU(i,j,k)=(grad(i,j)+grad(j,i))
!            TotalBulkViscosity = 0d0
        else
            TotalBulkViscosity = 0d0
        end if
    
    end function TotalBulkViscosity

end module DataPerfectGas